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ABSTRACT 

From pulsar scintillations we infer the presence of sheet-like structures in the ISM; it has been 
suggested that these are current sheets. Current sheets probably play an important role in heat¬ 
ing the solar corona, and there is evidence for their presence in the solar wind. Such magnetic 
discontinuities have been found in numerical simulations with particular boundary conditions, 
as well as in simulations using an incompressible equation of state. Here, I investigate their 
formation under more general circumstances by means of topological considerations as well 
as numerical simulations of the relaxation of an arbitrary smoothly-varying magnetic field. 
The simulations are performed with a variety of parameters and boundary conditions: in low, 
high and of-order-unity plasma-y^ regimes, with periodic and fixed boundaries, with and with¬ 
out a friction force, at various resolutions and with various diffusivities. Current sheets form, 
over a dynamical timescale, under all conditions explored. At higher resolution they are thin¬ 
ner, and there is a greater number of weaker current sheets. The magnetic field eventually 
relaxes into a smooth minimum energy state, the energy of which depends on the magnetic 
helicity, as well as on the nature of the boundaries. 

Key words: {magnetohydrodynamics) MHD - pulsars: general - turbulence - ISM: kinemat¬ 
ics and dynamics - ISM: structure - ISM: magnetic fields 


1 INTRODUCTION 

In this introduction pulsar scintillations, coronal heating, and some 
other relevant astrophysical contexts are reviewed. Afterwards in 
sections 2 and 3 I describe how we might expect magnetic discon¬ 
tinuities to form and present some simulations demonstrating this, 
before summarising in section 4. 

1.1 Pulsar scintillations 

Pulsar scintillations, a time variability observed in essentially all 
pulsars, were discovered shortly after the discovery of pulsars 
themselves (Scheuer 1968). They are produced by variations in the 
refractive index of the interstellar medium (ISM), which in turn 
depends on the density of free electrons rie. Lithwick & Goldre- 
ich (2001) proposed that the density variations are the result of 
Kolomogorov turbulence, but this picture is apparently inconsistent 
with reeent observations. Stinebring et al. (2001) found parabolic 
shapes in dynamic spectra when viewed in Fourier space, which are 
best interpreted as being produced by a small number of scattering 
screens in the line of sight (Walker et al. 2004; Hill et al. 2005; 
Cordes et al. 2006; Brisken et al. 2010). 

To produce the observed time delays we need path-length dif¬ 
ferences corresponding to seattering angles of order 1-100 mas. 
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The scattering angle is proportional to the density contrast Ane. 
If the scattering comes from sheets aligned at arbitrary angles to 
the line of sight, then we require a contrast Ane ~ lO^ne, which 
seems difficult to produce by any obvious physical process. Fortu¬ 
nately though, Snell’s law tells us that if the angle of incidence is 
small (i.e. grazing incidence) then the scattering angle is inversely 
proportional to the angle of incidence. Consequently we ean have 
Ane ~ ne as long as the angles of incidenee are of order 10“^. 
This is much easier to produce with plausible physical processes. 

Given that any mechanism to produce sheet-like structures is 
likely to produce them throughout the ISM, but that observations 
require only one or a few scattering structures at any one time, 
this grazing incidence picture has the advantage that if the required 
grazing angle is 10“^, only a fraction 10of sheets in the line 
of sight will produce appreciable scattering. To best explain obser¬ 
vations, they should be of order 0.1 AU in thickness, separated by 
perhaps 0.1 pe. An alternative geometry of filaments, as seen in 
cosmological dark matter simulations, not only requires alignment 
in two directions rather than one direction, but it is also physically 
difficult to explain, as self-gravitating objects for instance. 

Goldreich & Sridhar (2006) proposed that these structures are 
current sheets. Pen & Levin (2014) constructed a model where the 
seintillations are produeed by density variations associated with 
‘ripples’ - magnetie waves propagating along current sheets. This 
can reproduce the shapes seen in the secondary dynamic spectra. 


© 0000 RAS 



2 


Jonathan Braithwaite 


1.2 The solar corona 

Current sheets are thought to play a major role in coronal heating 
(Parker 1972, 2012 & refs, therein; Galsgaard & Nordlund 1996). 
As relevant to the topic of this paper, there are two differences be¬ 
tween the corona and the ISM: while the ISM has a plasma-of 
order unity the solar corona has a low-^^ and its magnetic field is 
subject to the force-free condition V x B aB. The other dif¬ 
ference is the boundary: motion in the corona is driven by convec¬ 
tive motion at the boundary, i.e. the photosphere. This motion has 
a timescale significantly longer than the dynamic timescale in the 
corona, so the coronal magnetic field evolves in quasi-static equi¬ 
librium. The ISM, on the other hand, lacks clearly defined bound¬ 
aries. 

In the ‘magnetic carpet’ picture of coronal heating, convective 
motion at the photosphere moves the footpoints of coronal mag¬ 
netic field arches around in uncorrelation, random directions. With¬ 
out topological reconnection, field lines become tangled around 
each other, with a complexity increasing with the passage of time. 
Parker (1972) argued that there is generally no new smooth force- 
free equilibrium solution for the magnetic field to relax into, so the 
new equilibrium must contain tangential discontinuities, which are 
presumably sites of reconnection and heating. Although this rea¬ 
soning does not apply close to a boundary with imposed smooth 
conditions, it should become asymptotically valid at large distance 
from the boundary. At a large distance from the boundary we are 
in a similar situation to the ISM - this argument could conceivably 
also apply to fluids without meaningful boundaries. 

A topic of debate is whether, when starting with a potential 
(i.e. current free) coronal field, it is necessary to move the foot- 
points a certain distance before discontinuities appear. If disconti¬ 
nuities appear as soon as the field is moved even slightly away from 
a potential configuration, and the discontinuities allow the field to 
relax efficiently back into a potential state, then the coronal field 
would be in a quasi-statically evolving near-potential state. This 
would be a serious problem for coronal heating, because moving 
the footpoints of a potential field requires no work, and no energy 
could be transferred from the convective motion into the corona. If 
however the discontinuities do not form immediately, or if the dis¬ 
continuities do not immediately produce rapid reconnection, then 
the coronal field would instead be in a quasi-static non-potential 
force-free equilibrium, presumably always slightly over the thresh¬ 
old for the appearance of current sheets (unless some instability 
threshold is crossed, as happens in fiares), and the moving foot¬ 
points would experience a magnetic drag force. 

Galsgaard & Nordlund (1996) performed simulations in which 
motion at the photosphere results in the formation of discontinuities 
in the corona, at which reconnection takes place. If the boundary 
driving is halted, the reconnection apparently stops but the discon¬ 
tinuities remain; reconnection at the discontinuities resumes once 
the driving is switched back on again. It would therefore be possi¬ 
ble for an equilibrium to contain discontinuities, but perhaps only 
with certain boundary conditions. 

Various authors have disputed this model. For instance Craig 
& Sneyd (2005) find in numerical experiments with boundary dis¬ 
placements that smooth equilibria can form. However, Low (2013) 
points out a boundary problem with their method of frictional re¬ 
laxation, leading to false solutions; damping via fiuid viscosity may 
be more physical than a friction force. 

^ The plasma-/3 is defined as the ratio of gas pressure to magnetic pressure: 
(3 = S7tPIB‘^. 


1.3 Other contexts 

Massive stars may also have magnetically active coronae. Stars 
above about 7 Mq have a convective layer just below the surface 
resulting from an iron-ionisation heat-capacity bump (see Cantiello 
et al. 2009). This layer could host a dynamo, and the field generated 
should easily reach the surface via buoyancy (Cantiello & Braith¬ 
waite 2011). These stars emit X rays, but with current observational 
data it is not possible to ascertain whether the source of X rays is 
close to the stellar surface, in which case there is probably solar- 
type magnetic heating, or further away, in which case it could be 
shock heating associated with the line deshadowing instability in 
the wind. However, various kinds of variability do point to activ¬ 
ity close to or on the stellar surface (Michaux et al. 2014 and refs, 
therein). 

There is evidence of discontinuities in the solar wind, inter¬ 
preted by Borovsky (2008) as flux tubes and by Li (2008) as cur¬ 
rent sheets. In an effort to explain the latter, Zhdankin et al. (2012, 
2013) performed simulations in reduced MHD, using an incom¬ 
pressible equation of state and a strong guide field. They find dis¬ 
continuities, and have some success in matching them statistically 
to observations. Greco et al. (2010) also find discontinuities in in¬ 
compressible simulations. Discontinuities were apparently already 
found in similar simulations by Muller & Biskamp (2000); unfortu¬ 
nately they are impossible for the reader to see because the results 
are presented only Fourier space. 

The Voyager probes recently fiew through the reverse shock 
of the solar wind (termination shock), and it appears that Voy¬ 
ager 1 may now be fiying through the contact discontinuity (he¬ 
liopause) into the (perhaps shocked) ISM (e.g. Schwadron & Mc- 
Comas 2013). Both probes have fiown through numerous current 
sheets. They may be the manifestation of the rotation of the Sun 
combined with the Sun’s magnetic field being tilted with respect to 
the rotation axis, but the irregular intervals between discontinuities 
speak instead in favour of magnetic reconnection events (Opher 
etal. 2011). 

Underdense X-ray bubbles in galaxy clusters, inflated by mag¬ 
netised AGN jets, often emit sychrotron radiation from cosmic rays. 
Often a considerable distance from the central AGN, the cosmic 
rays must presumably be produced in situ, perhaps in magnetic re¬ 
connection sites. The presence of these would not be surprising: as 
Braithwaite (2010) and Gourgouliatos et al. (2010) have argued, the 
magnetic field in a bubble should undergo plenty of reconnection 
over a long timescale. 


2 FORMATION OF DISCONTINUITIES IN MHD 

In this section the formation of current sheets is examined from 
the point of view of magnetic topology, using the ISM as an ex¬ 
ample context. The ISM is stirred amongst other things by super¬ 
novae, by stellar winds and ionising radiation, by the galactic disc’s 
self-gravity, shear fiow and associated instabilities. These sources 
drive motions on a variety of scales, and most of them are transient; 
in addition, the dynamic (acoustic or Alfvenic) timescale is often 
smaller at the stirring length scale than the intervals between stir¬ 
ring events. It is natural to expect therefore that at least locally the 
ISM relaxes towards some equilibrium state which evolves quasi- 
statically in response to changing external conditions. This state 
would persist until it is destroyed by some event. 

Note that the formation of current sheets was apparently men¬ 
tioned by Biskamp (1997), as part of a discussion of turbulent cas¬ 
cades. 
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2.1 Topological relaxation 

Locally one can consider what should happen after a stirring event. 
Initially, the terms present in the momentum equation are unbal¬ 
anced. The momentum equation and the sizes of the various terms 
are: 

^ = --VP + ^(VxB)xB + ( 1 ) 

at p 4:7Tp 

jt2 2 2 

U —C 3 va -j- 

p L 

where u, P, p, B and v are the velocity, pressure, density, mag¬ 
netic field and viscosity, where U and L are typieal fiow speeds 
and length scales, Ap is a typical density variation, and Cg and 
VA are typical sound and Alfven speeds. To get the second line, 
I have multiplied by L and assumed that U ^ L/T where T is 
a typieal timescale. In the ISM we have approximate equipartition 
between kinetic, thermal and magnetic energies, and the first three 
terms in the momentum equation are all eomparable in size. The 
last term, whose size vUjL can be written as /Re where Re is 
the Reynolds number, is enormously smaller than the other terms 
in almost all astrophysical contexts (one exeeption being the gas in 
galaxy clusters; see e.g. Schekochihin & Cowley 2006). This leads 
inevitably to nonlinear dynamics. 

In the absence of driving the system will head towards an en¬ 
ergy minimum. Kinetie and magnetie energy is dissipated into ther¬ 
mal energy via shocks and reconnection as well as by diffusion, at 
least on the smaller scales. As can be seen by comparing the sizes 
of terms in (1), this will happen on the dynamic timescale. The end- 
state is then an equilibrium where the pressure gradient and Lorentz 
forces balance each other: 

VP=^(VxB)xB. (2) 

47r 

The Lorentz force is perpendicular to the magnetie field, and in 
equilibrium must also be perpendicular to the isobars, meaning that 
field lines reside on surfaces of constant pressure. Picking an arbi¬ 
trary starting point and following a field line, one visits every point 
on a surfaee of eonstant pressure - one speaks of a surface-filling 
magnetic field. A general vector field, however, is volume-filling, 
meaning that following a field line one eventually visits every point 
in spaee. Generally, one would expeet that when a stirring event 
has finished, the ISM will be left with a volume-filling magnetie 
field^. In almost all astrophysieal eontexts the Ohmic timescale is 
very long, and fiux-freezing should hold to a high degree. How then 
is it possible to reach the rather special surface-filling topology 
required by equilibrium? The answer presumably involves small 
length scales, on which magnetic diffusivity can have a signifieant 
effect on a dynamic timescale, resulting in topological reconnec¬ 
tion. In other words, spontaneous formation of diseontinuities, or 
current sheets. Their thiekness presumably depends on magnetic 
diffusivity, but this need not affect the global behaviour, just as the 
viscosity of a gas affects the thickness of a shock but has no effect 
on conditions on either side. 

Note that the situation is different in two dimensions, where 
there is no topological barrier to equilibrium: in order to satisfy 
V • B every field line necessarily joins up with itself, and all that 
is required to reach an equilibrium is to equalise pressure along 


^ As long, of course, as the stirring also changes the previously established 
topology, which there is no reason to dispute since one expects reconnection 
during stirring as well as during relaxation. 


each line. Interestingly, Gruzinov (2009) (using periodie bound¬ 
aries) demonstrated that this equilibrium generally eontains dis¬ 
continuities, even though the initial conditions are smooth. In three 
dimensions, ArnoTd (1986) suggests that in the case of infinite eon- 
duetivity, the equilibrium should consist entirely of discontinuities. 
Obviously with high but finite conductivity something else hap¬ 
pens: perhaps a slowly-evolving equilibrium forms, containing dis¬ 
continuities separated by smoothly-varying regions, or perhaps an 
equilibrium forms with no discontinuities at all. The outcome could 
well depend on the boundary conditions. 

2.2 Magnetic helicity 

Magnetic helicity is a global quantity defined as the volume integral 
of the scalar product of the magnetie field with its vector potential: 
H = {1/ Stt) f A'B dV. In the case of infinite conductivity, the he¬ 
licity of a volume is conserved as long as no magnetic fiux passes 
through the boundaries (Woltjer 1958). It a topological property 
of the field which cannot be changed by moving the fiuid around 
while the fiux is frozen into it. It has units of energy x length and 
so is present (in a hand-wavy sense) more in the larger struetures 
than is the energy - and it seems, at least in many situations, to 
be approximately conserved when there is a small but finite mag¬ 
netic diffusivity, even whilst magnetic energy is being dissipated on 
small scales. This property has been demonstrated in many eontexts 
from the laboratory (Hsu & Bellan 2002) to the solar corona (Zhang 
& Low 2003). In section 3.5 it can be seen that helieity conserva¬ 
tion improves with deereasing magnetie diffusivity. Consideration 
of dimensions gives us 

\H\ = IE (3) 

where I is the some length scale and E is the magnetic energy. In a 
simple eonfiguration, I should be comparable to the size of system. 
For a given value of helicity the energy is lowest for a greater I 
and so equilibria should preferentially form with a high I and a 
simple, large-scale structure. However there may in principle be 
local energy minima with smaller I and greater E. 


3 NUMERICAL SIMULATIONS 

We are interested in how an arbitrary smoothly-varying magnetic 
field evolves in the absence of external driving and relaxes towards 
equilibrium. The setup here is a cube of side L eontaining gas with 
an ideal-gas equation of state and ratio of heat capacities 5/3. Ex¬ 
cept for the simulations described in section 3.3, boundaries are pe¬ 
riodie. In the initial conditions, we have uniform pressure and den¬ 
sity, zero velocity, and a magnetic field produced from a Fourier 
transformation of random phases in three-dimensional wavenum¬ 
ber spaee, at a range of wavenumbers from k = 27t/L (i.e. the 
smallest possible wavenumber) to Qtv/L, with a k~^ power law. 
The detail of this is not important; what is important is that this is 
a volume-filling magnetic field, which will require topological re¬ 
connection to reach a surface-filling equilibrium, and that the field 
contains only relatively large length scales. 

The code used is the STAGGER Code (Nordlund & Galsgaard 
1995; Gudiksen & Nordlund 2005), a Cartesian finite-difference 
MHD code with fifth- and sixth-order spatial derivatives and in¬ 
terpolations and third-order time-stepping. It has a hyperdiffusive 
scheme whieh damps structure near the spatial Nyquist frequeney 
whilst preserving structure on larger scales. The hyperdiffusion is 
switehed on for most of the simulations, but switehed off for a few. 
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where a slightly more quantitative examination of the effect of dif¬ 
fusion is requried. 

In the next section I describe a set of simulations differing 
only in their random initialisation. In subsequent sections I describe 
additional simulations with various different parameters. 

3.1 Plasma-/3 of order unity 

In this section, I present results from simulations run at a resolution 
of 128^ with hyperdiffusion switched on. The magnetic Prandtl 
number Prm = y/r] = 10 (where v and 77 are the viscosity and 
magnetic diffusivity), the reason being that the ISM, solar corona 
and essentially all other low-density astrophysical environments 
have high Prm and 10 is computationally easily achievable. The 
initial strength of the field is set such that the plasma-has a mean 
value of 1 / 2 , which, as we shall see below, results in an equilibrium 
with [3 ~ 1, as in the ISM. Several simulations were performed 
with these parameters, with different random seeds used to create 
the initial conditions. Note that the time units used in this paper are 
normalised such that an Alfven wave takes on average 1/2 time 
units to travel a distance L/27r in the initial conditions; this time 
becomes somewhat longer as the field decays. In other words, the 
time unit is roughly equal to an Alfven timescale. 

Fig. 1 shows snapshots from one of these simulations at three 
points in time, and fig. 2 shows snapshots from three more simu¬ 
lations. Multiple current sheets form in each case; the thickness of 
these sheets is just a few grid spacings. Following on from the ideas 
mentioned in section 2 . 2 , it is informative to look at the evolution 
of magnetic energy and helicity in the simulations; this is shown 
in fig. 3. The simulations begin with equal energies but different 
helicities. Clearly, helicity is approximately conserved, and energy 
drops to a level given by 

E = (4) 

where /cmin = 2t^ j Lis the minimum wavenumber in a box of size 
L. Note that one simulation (marked 1 in fig. 3) has zero initial 
helicity; its initial conditions are created by adding together two 
other magnetic fields in the correct ratio. This simulation displays 
the most reconnection activity, and the simulation with the greatest 
helicity (marked 4) displays the least. This is not surprising as the 
lower the helicity, the more magnetic energy is dissipated. Note that 
because of the ‘vanishing force-free field theorem’ the equilibrium 
end state is not force-free. (For the 3-line proof of this theorem see 
Roberts 1967, also reproduced in Spruit 2013.) 

What is actually happening? At first, since the Lorentz force 
is not balanced by anything, magnetic energy is converted into ki¬ 
netic. This happens on the Alfven timescale, as we see in the up¬ 
per panel of fig. 3. There is some sloshing around between mag¬ 
netic, kinetic and thermal energy - note how the magnetic energy 
rebounds a little after its initial drop. The average sonic Mach 
number, starting from 0 , will peak at somewhat less than unity 
at this stage. Tangential discontinuities appear: surfaces parallel 
to the magnetic field, with a field of different direction and mag¬ 
nitude on either side. Between these discontinuities, the field is 
smoothly varying. Over a dynamical timescale the number of cur¬ 
rent sheets increases until the situation becomes surprisingly com¬ 
plex: although only large length-scales are present in the initial con¬ 
ditions, eventually there are only relatively small distances between 
neighbouring current sheets. Current-sheet activity peaks and there¬ 
after the number and strength of current sheets falls as an equilib¬ 
rium is approached asymptotically. Note that current sheets appear 
even in simulation 4, which begins with a helicity equal to 90% 



time 



time 

Figure 3. The magnetic energy (thick lines) and helicity (thin lines) against 
time in four simulations, which differ only in their random initialisations, 
each plotted with a different line style. The upper panel follows the initial 
evolution, and the lower panel goes to later times. Helicity has units of 
energy times length so actually plotted on these axes is helicity divided by 
the length scale L/27r. Note how helicity is approximately conserved, and 
how the energy drops from its initial value of 498 towards the level set by 
the initial helicity. These simulations are all run at a resolution of 128 with 
Prm = 10 and hyperdiffusion. Note that simulation 4 begins with zero 
helicity. 


of the maximum, and which is therefore already rather close to an 
equilibrium. This may be relevant to the debate regarding whether 
discontinuities only form once the configuration is some distance 
from its energy minimum. A more thorough investigation of this 
point is left for the future. 

The discontinuities allow reconnection to take place on the 
dynamic timescale, or rather, somewhat more slowly: in fig. 3 we 
can see that after the initial ‘sloshing’ phase and once the current 
sheets have formed, the evolution slows down by around an order 
of magnitude. Such a factor of 10 between the Alfven speed and the 
reconnection speed is fairly standard in the literature (e.g. Eisner & 
Lamb 1984; Ikhsanov 2001). 

Note that at the beginning of these simulations, since the 
plasma-is of order unity, there is some sonic motion and weak 
shocks are present. Consequently, some discontinuities appear of 
the kind which have magnetic flux passing through them. The ma¬ 
jority of discontinuities are however of the tangential kind, and after 
one or two dynamical timescales have passed, all discontinuities are 
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Figure 1. Snapshots at t = 1.2, 4.3 and 18.3, from left to right, of the simulation marked 2 in fig. 3. The volume rendering represents current density, with the 
same thresholds in all plots: green, blue and purple represent medium, high and very high current density. Current sheets form on the dynamic timescale. The 
second snapshot is at the time of maximum activity, after which the number and intensity of current sheets drops as equilibrium is approached. The simulations 
was run at resolution 128^ with Prm = 10 and hyperdiffusion, and began with a mean plasma-/3 of 1/2. Note that the mottling effect on the current sheets 
comes from the visualisation procedure. 



Figure 2. Snapshots of three simulations identical to that in fig. 1 except for the random seed used to generate the initial conditions; they are marked 1, 3 and 
4 in fig. 3. The snapshots are taken in each case approximately at the time of maximum activity: att = 2.0, 3.3 and 4.2, from left to right; the magnetic 
helicity in the initial conditions increases in that order. The simulation on the left has zero initial helicity. Note the rolled features on the right-hand-side of that 
snapshot - these concentric current sheets annihilate each other shortly after appearing. 


tangential. In the following section a simulation is presented with a 
high plasma-/^, from which shocks are obviously absent. 

3.2 The high- and low-/3 regimes 

The simulations in the previous section had a plasma-/3 of approx¬ 
imately unity, as in the ISM. It would now be interesting to see 
whether the same behaviour exists in the regimes with high and 
low plasma-/3; previous results using an incompressible equation 
of state suggest that current sheets should form in the high-(3 case. 
To this end, simulations are run with the initial plasma-/3 set to 15 
and 1 /60, which evolve to around 30 and 1/30. Ohmic and viscous 
heating are switched off for the latter, in order to remain in the low- 
/3 regime^. Fig. 4 shows snapshots from simulations at high and 
low P, using the random initialisation corresponding to simulation 
2 in section 3.1. 

The magnetic fields in these simulations have taken different 
forms, but are qualitatively similar: it is probably safe to conclude 

^ Clearly this is unphysical; real astrophysical plasmas can only remain in 
the low-/3 regime with the help of physics not present in these simulations, 
such as radiative cooling, gravity and particular boundaries. 


that current sheets form in all three regimes. In principle the be¬ 
haviour might be different at values of /3 more than a factor of 30 
away from unity, but there is no particular theoretical reason to ex¬ 
pect this. 

3.3 Boundary conditions 

The role of boundary conditions in the formation of discontinuities 
is a topic of some debate, particularly in the context of the solar 
corona, into which magnetic energy is continually injected via mo¬ 
tion at the boundary, i.e. the photosphere. It is therefore a useful 
exercise to repeat the simulations described above with different 
boundary conditions. Those most akin to the coronal context would 
be a fixed boundary in one of the three dimensions, as envisaged in 
Parker’s original model. 

In fig. 5 we see snapshots from three simulations, all of which 
use the same initial conditions as simulation 2 from section 3.1. 
One has a fixed boundary in the x direction, one has fixed bound¬ 
aries in all directions, and one has a ‘pseudo-vacuum’ boundary 
in one direction. The latter condition has symmetric p, P, and 
B±, antisymmetric u± and By, and is essentially a simple way 
to get similar results to a potential field extrapolation. In addition, 
simulations were run with each of these boundary conditions in 
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Figure 4. Snapshots of two simulations at f 0.6 and 14.5, with identical initial conditions except for the mean plasma-/?, which is 1/60 and 15 (left and 
right). The volume rendering represents current density. Clearly, the spontaneous formation of discontinuities is a phenomenon common to both high- and 
low-/? regimes. In both simulations, the spatial resolution is 128^, with Prm = 10 and hyperdiffusivity switched on. 


the high- and low-/? cases (the same initial conditions as in sec¬ 
tion 3.2). Current sheets always form, but are in different places 
in each case; as with periodic boundaries, they appear over a dy¬ 
namical timescale. The current sheets then become weaker as an 
equilibrium is approached asymptotically. In the simulations with 
periodic boundaries, it seems that the equilibrium is smooth (cor¬ 
responding to minimum energy), but with fixed boundaries it is not 
completely clear whether the equilibrium contains discontinuities 
too weak and/or numerous to be seen in the simulations. 

In fig. 6 are plotted the energy and helicity in the three simula¬ 
tions with initial /? = 1/2, together with simulation 2 from section 
3.1, which has periodic boundaries. Fixed boundaries significantly 
slow the drop in magnetic energy relative to the case with periodic 
boundaries, and the end state is different, with a higher energy. 
Pseudo-vacuum boundaries do not conserve helicity, as expected. 
A more detailed study of the role of boundaries is described in a 
forthcoming paper. 

3.4 Frictional relaxation 

In many studies of relaxation into a minimum energy state, the fiuid 
is not free to slosh around, but is instead forced to move directly 
towards an equilibrium. This is often done by assigning the fiuid 
velocity as the sum of the pressure gradient and Lorentz forces 
times some multiplier. In other words, the momentum equation (1) 
is modified: the left hand side changes from du/dt to u/rfric where 
Tfric is a constant timescale, and the viscous term is dropped. This 
simplifies the relaxation process; the system finds an equilibrium 
much faster. However, according to Low (2013), this method can 
lead to the wrong equilibrium being reached. 

To change fundamentally the nature of the momentum equa¬ 
tion within the framework of the MHD code used here would be 
time consuming, but fortunately there is a simpler solution to mimic 
this: the addition of a friction force —u/rfric on the right-hand 
side of the momentum equation (1). We have something like criti¬ 
cal damping if Tfric is comparable to the dynamic timescale, or in 
other words if the friction term on the RHS and the inertia term 



time 


Figure 6. Energy (thick lines) and helicity (thin lines) in four simula¬ 
tions with identical initial conditions but different boundaries. The pseudo¬ 
vacuum boundary allows helicity to decay on the dynamical timescale, as 
expected. The two simulations with fixed boundaries display an approach 
to equilibrium at constant helicity, but this happens rather more slowly than 
in the case with periodic boundaries in all three directions. All simulations 
were run at a resolution of 128 with Prm = 10 and hyperdiffusion, with an 
initial average plasma-/? of 1/2. 

on the LHS of the momentum equation (1) are comparable. If Tfric 
is smaller than the dynamic timescale, the fiuid is ‘overdamped’ 
and inertia becomes irrelevant; this is the regime closest to the 
pure frictional relaxation, we which wish to investigate. Calling 
the timescale on which the field would otherwise evolve to, which 
as we saw above is comparable to ta or a little longer, then in 
the regime Tfric ^ tq we expect the evolution to take place on a 
timescale tq jTfric - This comes from assuming balance between the 
Lorentz force and the friction force, as opposed to balance between 
the Lorentz force and inertia as we had previously. 

Simulations were run in the high-, low- and of-order-unity-/? 
regimes (as described in section 3.2) with friction force timescales 
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Figure 5. Snapshots from three simulations using the same initial conditions and parameters as simulation 2 from section 3.1, but having different boundary 
conditions. Left: fixed in one dimension (left and right sides, i.e. dXx = -^Ll2, as plotted here). Centre: fixed in all three dimensions. Right: pseudo-vacuum 
in one dimension (also left and right sides). The snapshots are at times t = 3.3, 3.6 and 2.5. Current sheets form in all cases, just as with periodic boundaries. 
All three are run at a resolution of 128 with Prm = 10 and hyperdiffusion, with an average initial plasma-/? of 1/2. 



Figure 7. A snapshot at time t = 60 from a simulation with a friction force 
with a timescale Tfric = 0.1, which is about a tenth of the Alfven timescale. 
Otherwise, this simulation is the same as the simulation 3 from section 3.1. 
In comparison to that simulation we have here fewer, and somewhat thicker, 
current sheets. 


in each case approximately ten times smaller than the Alfven 
timescale. In fig. 7 we see a simulation otherwise identical to simu¬ 
lation 3 from section 3.1, so with initial plasma-/? set to 1/2. There 
are fewer current sheets, which is presumably because the slower 
relaxation towards equilibrium means that some current sheets have 
already come and gone before others form, rather than all form¬ 
ing at the same time. Also, the reduction in reconnection speed 
increases the thickness of the current sheets, which is not at all 
surprising. 

We see essentially the same qualitative behaviour: the forma¬ 
tion of current sheets on a dynamic timescale, in all three plasma-/? 
regimes. Inertia plays no role in these simulations, so we see that 
the formation of current sheets is essentially a matter of topology. 
Importantly though, the addition of a friction force does result in a 
quantitatively completely different end-state, albeit with the same 
magnetic energy. In some contexts this difference might be quite 
important. 



Figure 9. Snapshot of the resolution 256 simulation at the same point in 
time as in fig. 8 (right-hand panels). The colour thresholds for current den¬ 
sity are slightly different from that figure, and field lines are shown in red. 


3.5 The effects of diffusivity and resolution 

A degeneracy is expected here since increasing the resolution is 
conceptually the same as decreasing the diffusivity. Or to be more 
precise, increasing the resolution allows a reduction in diffusivity. 
The simulations presented in previous sections were all run at a res¬ 
olution of 128^, with Prm = 10 and hyperdiffusivity switched on, 
with the minimum magnetic diffusivity required to reliably remove 
unpleasant zig-zags. 

First, we can look at simulations with resolutions of 64^ and 
256^, i.e. with higher and lower effective diffusivities, but which 
are otherwise identical. Fig. 8 shows snapshots from simulations at 
all three resolutions, using the random initialisation corresponding 
to simulation 2 from section 3.1. Clearly, the strongest features are 
visible in all three frames. At higher resolution, the current sheets 
are somewhat thinner, and there is an abundance of weaker features 
which are not present at low resolution. Fig. 9 shows the 256 reso¬ 
lution run (same snapshot) with field lines. 

It is informative to look at the distribution of current density 


© 0000 RAS, MNRAS 000, 1-10 



















8 Jonathan Braithwaite 



Figure 8. Snapshots of three simulations att ^ 2.5, with identical initial conditions but with spatial resolutions of 64, 128 and 256 (left to right). Top row: 
the volume rendering represents current density. Bottom row: slices through the simulation with current density represented by a similar colour-coding as in 
the 3D plots: red is below-average, then increasing current density is green, blue, then purple. Whilst the larger features are present at all three resolutions, 
there are additional weaker current sheets at higher resolution. All simulations have an initial mean plasma-/3 of 1/2. 



Iogio J 


logio J 


log^Q J 


Figure 10. Histograms of current density in the three simulations shown in hg. 8: left to right, at resolutions 64, 128 and 256, each at various points in time. 
In these simulations Prm = 10. Note that over the dynamical timescale regions with very high current density develop. These then gradually fade away as 
the magnetic field approaches its minimum energy state, the current sheets becoming weaker. This effect is more pronounced at higher resolution, where the 
discontinuities are thinner and the current density correspondingly greater. 


in the simulations. In fig. 10 histograms of current density are plot¬ 
ted for the three simulations in fig. 8 at various points in time. We 
see that after only a short time, some modest fraction of the vol¬ 
ume hosts very high current density. This is the tangential discon¬ 
tinuities. After their formation, they fade in intensity as the field 
evolves towards its energy minimum - the magnetic field vectors 
on either side of a discontinuity gradually approach each other. 

One can also investigate the effect of varying the magnetic 
diffusivity, or in other words, varying the magnetic Prandtl number 
while holding the viscosity constant. Simulations were run other¬ 


wise identical to simulation 2 from section 3.1, which had Prm=10, 
with Prm =3 and 1. To also allow a more quantitative study, three 
further simuations were run with normal, rather than hyper- diffu¬ 
sion, with the same values of Prm - The effect is essentially the same 
as that of changing the resolution: lower diffusivity leads to thinner 
current sheets and more weaker sheets. In fig. 11 we have at the top 
a plot of energy and helicity in these six simulations, and below, the 
timescale of helicity decay, calculated as Thei = (d In |iT|/dt)“^. 

In these simulations with normal, not hyper- diffusion, helicity 
conservation begins to fail significantly. Helicity decays - it never 
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time 



time 


Figure 11. Upper panel: the magnetic energy (thick lines) and helicity (thin 
lines) against time in six simulations. All have the same initial conditions 
and are run at resolution 128^; the difference is in the diffusivity: three are 
run with hyperdiffusion and three with standard diffusion, each with mag¬ 
netic Prandtl numbers of 10, 3 and 1. The viscosity is the same in all simu¬ 
lations. Lower panel: For the same simulations (with the same line styles), 
the evolution timescale of the helicity defined as ruei = (d In \ H\/dt)~^. 
Clearly, the lower the magnetic diffusivity, the better the helicity conserva¬ 
tion. The hyperdiffusion scheme also improves helicity conservation. 


grows - and it decays faster in simulations with higher magnetic 
diffusivity. Furthermore, in the simulations with ‘normal’ diffusion, 
there is a visible quantitative correlation, namely that rate of helic¬ 
ity decay is proportional to magnetic diffusivity. With hyperdiffu¬ 
sion, helicity conservation is better and the same trend is present, 
albeit with a gentler scaling with diffusivity. This trend seems to 
hold at all times, right from the beginning when everything is mov¬ 
ing around on the Alfven timescale, through the peak of reconnec¬ 
tion activity to later times when the current sheets have become 
much weaker and fewer. At late times the helicity decay timescale, 
at least in the case with normal diffusion, approaches an asymptotic 
value equal simply to the global Ohmic timescale. 

Another effect of the magnetic Prandtl number is on how the 
magnetic energy is dissipated: we can see in fig. 12 that Ohmic 
dissipation accounts for a greater fraction of total dissipation when 
Prm is smaller. In astrophysical contexts with high Prm, such as the 
ISM where its value is around 10^^, that would mean that almost 
all of the energy is dissipated viscously. 



time 

Figure 12. For the same six runs as in fig. 11, the Ohmic cumulative en¬ 
ergy dissipation as a fraction of the sum of Ohmic and viscous cumulative 
energy dissipations. As expected, the larger one diffusivity with respect to 
the other, the greater energy it dissipates. There is also some bias towards 
Ohmic diffusion. 

4 CONCLUSIONS AND DISCUSSION 

Simple simulations were performed in a cubic computational box 
with an ideal-gas equation of state, uniform initial density and pres¬ 
sure, a smoothly-varying, large-scale initial magnetic field, and no 
driving of any kind. From these simulations we find the following: 

• Tangential discontinuities form on a dynamical (Alfven) 
timescale, and topological reconnection then proceeds on a 
timescale roughly a factor 10 greater than that; gradually the dis¬ 
continuities become fewer and weaker as an equilibrium is ap¬ 
proached. The eventual equilibrium is smooth, containing no dis¬ 
continuities. 

• Quantitatively, the behaviour is the same in all three plasma-^^ 

regimes: ^ ^ ^ 1 and ^ ^ 1. 

• Current sheets form in simulations with fixed and pseudo¬ 
vacuum boundaries, as well as periodic. 

• The addition of an artificial friction force makes no qualitative 
difference, demonstrating that inertia is irrelevant in the formation 
of current sheets. Importantly though, the end-state is quantitatively 
changed by the friction force. 

• At higher resolution (or lower diffusivity) the discontinuities 
are thinner. Whilst the strongest discontinuities are the same at high 
and low resolution, at high resolution there is a greater number of 
weaker discontinuities. 

• The magnetic helicity H is approximately conserved during 
the whole process, and at equilibrium, the magnetic energy E is 
given hy E = kminH where kmin = 27r/L is the minimum possi¬ 
ble wavenumber in a computational domain of side L. 

• In simulations with large helicity, fewer discontinuities form 
than when helicity is small. 

• Helicity conservation is better at higher resolution; with nor¬ 
mal diffusion (i.e. not hyperdiffusion) the rate of decay of helicity 
is proportional to magnetic diffusivity. 

Tangential discontinuities have already been found in stud¬ 
ies of the solar corona, where they are thought to be responsible 
for much of the heating. The focus of work in this context has 
been in forming discontinuities by jiggling the field lines around 
at the boundaries. Parker (1972) argued that there are generally no 
smooth solutions to such displacements at the boundaries, provided 
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that the displacements are of sufficient magnitude and complex¬ 
ity, a condition it is hard to imagine not being met in reality. All 
studies of the corona have assumed boundary conditions where the 
fiux is frozen into the boundary and subject to motion parallel to 
the boundary, mimicking the convective motion we see at the solar 
photosphere. Usually the field is assumed to be force-free, as the 
corona has a low plasma-^5. 

Tangential discontinuities have also been found in simulations 
with periodic boundaries and an incompressible equation of state in 
‘reduced MHD’, allowing perturbations to a strong uniform guide 
field (Zhdankin et al. 2013), and also without a guide field (Greco 
et al. 2010). 

From the simulations presented here, we see that the forma¬ 
tion of current sheets is very general: they form in both high- and 
low-y5 regimes, with various boundaries, with and without a large 
friction force, and from initial conditions which are topologically 
either rather close to the eventual equilibrium or very distant. In 
light of this, one should expect discontinuities to occur in many 
astrophysical contexts, including in the ISM, where they could ex¬ 
plain the observed properties of pulsar scintillations. 

In the future one might like to explore more thoroughly the ef¬ 
fect of the initial conditions, of which only one type was used here: 
large scale and smooth. It is conceivable that there are initial con¬ 
ditions which do not produce discontinuities as the magnetic field 
relaxes, but these might have to be either ‘special’ in some sense, 
or very close to the eventual equilibrium, a situation it is difficult to 
imagine arising in reality, except perhaps where an existing equi¬ 
librium is ‘stirred’ a just little and then allowed to relax back again. 
Another possibility is that there is some possible constellation of 
initial and boundary conditions which produces a equilibrium con¬ 
taining discontinuities. 

Finally, a comment on the observational effects in terms of 
pulsar scintillations. In the simulations, the gas pressure in the cur¬ 
rent sheets is generally somewhat different from the surroundings. 
Of interest for pulsar scintillations is the electron density rather 
than the pressure, about which these simulations tell us nothing 
directly, owing to their lack of the relevant thermal physics: they 
do contain Ohmic heating, but there is no radiative cooling, the 
thermal conduction is not particularly realistic, and the magnetic 
Prandtl number is far smaller than in the ISM. The current sheets 
are heated by reconnection, and will cool on some radiative cooling 
timescale which should be around 10^ yr at a density of 1 cm“^. 
The longevity of these current sheets is somewhat greater than the 
dynamical timescale; taking the sheets to be separated by distances 
of around 0.1 pc (inferred from observations) and taking the Alfven 
speed to be 10 km s“^, this is also around 10^ yr. Given that these 
timescales are comparable, it is not obvious whether the gas density 
in current sheets should be lower or higher than in the surroundings. 
For pulsar scintillations though the sign of the density difference is 
not important, only the magnitude, which will be a factor of or¬ 
der unity, or somewhat less than unity in locations where the ISM 
has been left relatively undisturbed for a long time and only weak 
discontinuities remain. 
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